function d_y = diff_t(t,y)
    n1 = length(t);
    d_y = [(y(2)-y(1))/(t(2)-t(1));(y(3:n1)-y(1:n1-2))./(t(3:n1)-t(1:n1-2)); ...
    (y(n1)-y(n1-1))/(t(n1)-t(n1-1))];
end